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Relaxation schemes for finding normal modes of nonlinear excitations are described, and 
applied to the vortex-spinwave scattering problem in classical two-dimensional easy-plane 
Heisenberg models. The schemes employ the square of an effective Hamiltonian to ensure 
positive eigenvalues, together with an evolution in time or a self-consistent Gauss-Seidel 
iteration producing diffusive relaxation. We find some of the lowest frequency spinwave 
modes in a circular system with a single vortex present. The method is used to describe 
the vortex-spinwave scattering (S-matrix and phase shifts) and other dynamical properties 
of vortices in ferromagnets and antiferromagnets, for systems larger than that solvable by 
diagonalization methods. The lowest frequency modes associated with translation of the 
vortex center are used to estimate the vortex mass. 

PACS numbers: 75.10.Hk; 75.40. Cx 



I. INTRODUCTION: VORTEX NORMAL MODES 

The spectrum of small-amplitude vibrations in the 
presence of a nonlinear excitation, such as a vortex, soli- 
ton, or other inhomogeneous magnetization, contains in- 
formation about the scattering properties, the transla- 
tional properties, and the internal modes of vibration and 
instabilities of that object. Any modes that acquire zero 
frequency, for example, as a parameter is changed, sig- 
nal a change of symmetry or phase transition, or perhaps 
a process such as magnetization reversal. Other modes 
may be coupled to the translation of the center of the ob- 
ject: such translation modes contain information about 
the effective mass of the excitation. If the modes are cal- 
culated by diagonalization methods, the size of the ma- 
trix and cpu time necessary to solve it very easily exceed 
the practical limits of any computer even for fairly small 
systems. Thus it is interesting to consider other methods 
for the calculation of at least the lowest frequency modes 
in the presence of an inhomogeneous state. 

In particular, classical magnetic models are known to 
support localized or at least partially localized nonlinear 
excitations, one of the most well-known examples being 
the vortices of two-dimensional (2D) models with easy- 
plane (XY) anisotropy: 

iJ = — — ^(•S'nS'i^+a + ^a^n+EL + '^'^'n'S'n+a) 
n.a 

-h-Y^Sr., (1) 
n 

where subscripts n and a label the lattice sites and dis- 
placements to the nearest neighbors, and < A < 1 
determines the anisotropy, and we have included an ap- 
plied magnetic field h. Quite generally, linearization of 



the associated equations of motion about a nonlinear ex- 
citation such as a vortex will lead to a normal mode or 
spinwave problem, with an operator M. producing the 
time evolution of a spinwave wavefunction ^': 

^i^ = M-^. (2) 

at 

The calculation of the spectrum of spinwaves in the pres- 
ence of a vortex by diagonalization was accomplished^ 
for small circular systems up to a radius R ~ 20a, where 
a is the lattice constant. The method was used to deter- 
mine vortices' spinwave spectrum for ferromagnets (FM, 
J > 0) and antiferromagnets (AFM, J < 0) on a lattice.^ 
In addition to usual scattering states, it was found that 
FM and AFM vortices possess an internal mode whose 
frequency goes to zero at a critical anisotropy value Ac, 
which is lattice dependent (« 0.7034 for square lattice). 
The mode is associated with the instability of out-of- 
plane vortices (those with 0) to change to in- 

plane vortices = everywhere) at strong easy-plane 
anisotropy (A < Ac).'^^^ Other modes found, which tend 
to have high intensity at the vortex core, are associated 
with the translation of the vortex center of mass, and 
can be related to a collective coordinate description for 
vortex dynamics. ^'^ 

Analytic calculations of the vortex normal modes, us- 
ing a continuum limit, have partially described the scat- 
tering spectrum. In the earliest calculations, the 
phase shifts of spinwaves scattered by FM in-plane vor- 
tices were calculated, using a Born approximation, for 
the exchange anisotropy model (1) considered here^ and 
for a similar FM model with site-anisotropy.^ The weak 
point of these calculations is that they rely on the choice 
of a short distance cutoff for the Born approximation in- 
tegrals. For the AFM in-plane vortices, Pereira et a/.^° 
were able to make an exact calculation of the scattering 
phase shifts, for the optical spinwave branch, without the 
need for a cutoff. In the continuum limit equations, AFM 
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vortices do not scatter the acoustic branch spinwaves; the 
associated phase shifts vanish. 

The vahdity of the continuum hmit and other approx- 
imations, however, can only be tested by comparison 
with numerical calculations. For example, Ivanov et al}^ 
found by analytic means that the vortices of the contin- 
uum AFM model for A near 1 should possess a true in- 
ternal local mode. It was shown only through numerical 
diagonalization that this mode is identical to the vortex 
instability mode mentioned above. Ivanov et al. were 
also able to calculate the scattering (S-matrix) by a nu- 
merical shooting solution for the continuum theory equa- 
tions of motion. The S-matrix was not calculated from 
the numerical diagonalization data, however, because the 
system size that could be solved was too small. 

For Ac < A « 1, the nonzero component of the out- 
of-plane vortex leads to a smooth structure decaying over 
the length scale 

r. = (3) 

with no singularity at the vortex center. Thus, the con- 
tinuum description used in Ref. 12 is quite accurate. At 
the other limit, A < Ac, the in-plane vortex spin field is 
singular at the vortex center. In this case, the contin- 
uum limit cannot describe very well the spin dynamical 
motions near the vortex core, without the introduction 
of some kind of lower radius cutoff. The length scale 
in the problem approaches zero, and derivatives of the 
static vortex spin structure on a lattice arc not small in 
the core region, violating the usual continuum limit as- 
sumptions about small spatial derivatives across one lat- 
tice constant. Nevertheless, even for the in-plane AFM 
vortices, a true local mode of the vortex is present, ^"^ 
and can be approximately described in continuum the- 
ory, but only if a cutoff is introduced. Therefore, for 
both FM and AFM in-plane vortices, it is interesting to 
consider the calculation of the S-matrix (or equivalently, 
phase shifts) for the lattice system, and compare with 
continuum scattering theory. 

It has been realized-'^^ that the translational modes that 
come from the spinwave spectra have unusual properties 
and c;an be used to calculate a vortex effective mass. Al- 
though only a limited size of system could be studied 
{R < 20a) it was found that the FM vortex mass in- 
creases faster than In R, for in-plane as well as out-of- 
plane vortices. AFM in-plane vortices also were found to 
have a mass increasing faster than Ini?, however, it ap- 
peared that the mass of AFM out-of-plane vortices may 
actually reach a finite limit for large R. It is also possi- 
ble that AFM vortex mass tends to zero for A — > 1, the 
isotropic hmit. 

More recently, a newer theory for the FM out-of-plane 
vortex mass and related collective coordinate description 
has been developed by Mertens et al.^ There, a vortex 
mass approximately independent of system size was de- 
termined, both by including effects due to an image vor- 



tex outside the studied system, and also, by using a dy- 
namical equation of motion higher than second order in 
time. This newer theory appears to describe well the indi- 
vidual vortex dynamics by assigning a mass, gyrovector, 
and a new higher order gyro-tensor to the vortex, con- 
sidered as a particle-like object. A careful application of 
the theory was made by Ivanov et al,^ using spinwave 
translation mode data obtained from another numerical 
scheme for finding some of the lowest modes. In Ref. 
7, Wielandt's version of inverse iteration procedure was 
used: an operator {M. — u}I)~^ was applied to an initial 
randomly chosen vector, which after many iterations con- 
verges to an eigenvector of A4 , provided the frequency lo 
is updated appropriately. It is also a type of relaxation 
procedure similar to what we describe here, with great 
speed and memory advantages over the diagonalization 
method. This vortex mass theory, however, does not ap- 
ply directly to the AFM vortices, nor to the in-plane 
vortices for both FM and AFM models. 

In this paper we present and apply some numerical 
methods for finding the normal modes of a magnetic vor- 
tex on a lattice, which are applicable to systems larger 
than those solvable by diagonalization schemes. We use 
these methods to analyze the FM and AFM spinwave- 
vortex scattering problems for in-plane vortices. In addi- 
tion to the calculation of the S-matrix and phase shifts, 
we present results on the translational modes and make 
estimates of in-plane vortex mass. Analysis of out-of- 
plane vortices will be considered elsewhere. 

The paper is organized as follows. First, in Sec. II, a 
short overview of relaxation schemes is given. In Sec. Ill, 
we describe the magnetic spinwave equations of motion 
in the presence of a nonuniformly magnetized configura- 
tion. This is followed in Sec. IV by a description of the 
relaxation schemes using time evolution and Gauss-Seidel 
iteration. In Sec. VI the schemes are used to calculate the 
vortex-spinwave S-matrix, for a single vortex in a circular 
system. After a summary of the translational dynamics 
of magnetic vortices, in Sec. VII we present and analyze 
the results obtained for the FM and AFM in-plane vortex 
mass. 



II. RELAXATION METHODS FOR 
EIGENFUNCTIONS 

The basic idea of these relaxation schemes is the fol- 
lowing: The eigenfrequencies of a linear spinwave prob- 
lem, such as Eq. (2), always come in positive/negative 
pairs, corresponding essentially to creation and annihi- 
lation operators, respectively. Thus, its time evolution 
leads to oscillatory behavior, starting from an arbitrary 
initial condition. It is interesting, however, to consider 
how the evolution can be made diffusive, such that high 
{absolute magnitude) frequency modes decay away faster 
than low frequency modes. For example, a Schrodinger 
equation with the replacement it r, becomes a dif- 
fusion equation in imaginary time r (because it is a 
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parabolic equation). It is well known that its imaginary 
time evolution will lead to the lowest frequency wave- 
function, i.e., the ground state of the associated Hamilto- 
nian. In the equations of motion (2), such a replacement 
it ^ T, unfortunately, does not lead to a diffusion- like 
equation, due to the presence of the two fundamental 
solutions (creation/annihilation operators) behaving as 
giiwfct _^ Q^<^kT (because the spinwave equations are es- 
sentially hyperbolic). There is a fundamental solution 
growing with r for every mode A;, leading to numerical 
instability in imaginary time. However, if one uses the 
square of the spinwave Hamiltonian to evolve forward 
in real time, then both fundamental solutions associated 
with mode k decay as e~'^'=*. Evolution forward in time 
will lead to the lowest frequency cigcnmodc. By com- 
bining with Gram-Schmidt orthogonalization, it is then 
possible to iterate this procedure and generate a small set 
of the lowest frequency eigenmodes. The method offers a 
drastic savings in memory needed compared to numerical 
diagonalization, at the expense of rather long cpu time 
due to its diffusive behavior. Also note, it is not essen- 
tial that the evolution be performed in time. We have 
also found that the eigenmodes can be found even more 
quickly by using a Gauss-Seidel iteration of the spinwave 
equations, where the unknown frequency uj within those 
equations is determined self-consistently during the iter- 
ation. 



III. 2D EASY-PLANE MAGNETIC MODEL AND 
SPINWAVE EQUATIONS 



(cos 6^ cos (t)^ 
cos 6*0 sin< 
sin 



(6) 



Then one could obtain the spinwave effective equations 
of motion by assuming a perturbation to this structure 
in the form. 



(7a) 



(7b) 



where the equations of motion are now linearized in terms 
of (/3n and 'dn- 

In Rcf. 1 a slightly different but equivalent description 
was used, which we follow here. The unperturbed spins of 
the static vortex structure are considered to define local 
quantization axes z^, specifically. 



(8) 



Then the perturbation of this structure involves fluctua- 
tions orthogonal to the Zn-axis, along two new local axes 
i„ and y^- The in-axis is taken to be along the direc- 
tion defined by the cross product Xn = x 5n, which is 
within the original xy (easy) plane. Quantities x^. y^, 
and Zn here refer to the original coordinates of Hamilto- 
nian (1). Then yn = -Sn x makes the last member of 
the local orthogonal set. The perturbation of the static 
vortex structure can then be expressed in terms of its 
spin components along these new local axes: 



We review briefly the derivation of the linearized spin- 
wave problem for normal modes on a vortex, or other 
nonuniform state. The nonlinear equations of motion 
from (1) are simple: 



S„ X 



h+ J-^- 



(4) 



n'=n-t-a 

where J is a constant diagonal exchange matrix: 

1 
J = J I 1 
A 



(5) 



and the spins are considered as column vectors. A vor- 
tex is a particular static configuration {S^} that solves 
this equation with a zero time derivative. The structure 
of in-plane and out-of-plane vortices has been described 
elsewhere using continuum theory^^ and in the presence 
of a lattice. ^'^ Because the spin length is a conserved 
quantity, there are effectively only two degrees of free- 
dom per spin, and it is usually convenient to describe 
the static structure using an in-plane angle and an 
out-of-plane angle in planar spherical coordinates: 



(9) 



A short calculation shows that these are related to the 
angular perturbation coordinates by 



(10) 



The variables and i^n relate to purely in-plane spin 
motions, while S*^ and ?9n measiirc the out-of-easy-plane 
tilting, relative to the unperturbed vortex. 

In Refs. 1 and 2 it is shown that under the assumptions 
S'^ <C 1, S'l ^ 1, S*^ ~ S', this notation leads to a usual 
linear time evolution problem, which we write here in the 
form: 



It 



E 

m=n,n' 



Myy 



or equivalently in components. 



(11) 



(12) 



m— n,n' 



where a and /3 range over the set {i, y}, and the site m 
must either be equal to site n or one of its neighbors, 
n' = n -|- a. The matrix elements M"^ are: 
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J^XX 

nn 



Mil 



nn 



0, 



^^■^ n,n' 



{h^ cos + hy sin + h 



n'=n+a 



JS m°sin(^° 
M%, = -JS [mlmi, cos(<^o - + Xpip^,] 



(13a) 


n 

(13b) 

(13c) 
(13d) 
(13e) 
(13f) 



For simplification, wc use the definitions of in-plane pro- 
jection and out-of-planc' magnetization, 



Pr, 

ml 



Vl -(5^7^)2= cos C 
S^^/S- 



■ sm 



9°. 



(14a) 
(14b) 



The operator on the right hand side of the Uncar equa- 
tions of motion, Eq. (12), is not Hermitian. (For example, 
consider the case A = 0, with all ^° = 0.) This led to 
complications in the previous numerical diagonalization 
calculations. Its eigenvalues are pure imaginary, in com- 
plex conjugate pairs; this results in the usual spinwavo 
solutions oscillatory in time. However, this implies that 
its square (which need not be Hermitian) has pure real 
eigenvalues, a fact wc will exploit below in Sec. IV. 

In Refs. 1 and 2, the variables 5^ and S^ were con- 
sidered as quantum operators, in a semiclassical picture. 
Then an eigcnmodc labeled by index fc, has a creation 
operator B^. that is an appropriate linear combination, 
with complex expansion coefficients, wl^,iu^^, 



(15) 



such that its time evolution is simple, with eigenfre- 
quency w^: 



Bl = iwkBl 



(16) 



The corresponding annihilation operator is B^. = {Bj,)*, 
since the spin operators are real. The local spin fluctua- 
tions are given from the inverse relationship. 



s^ = insJ2{<r.Bk 



(17a) 



Si = -ins J2 {<r.Bk - wi*Bl) , (17b) 



Note that the italic superscripts on the w'^, a = 1 ,2 in 
this article should not be confused with powers, because 
this is a linear problem. 

The collection of coefficients describes the wavefunc- 
tion, and satisfies 



lUJk 



w 



E 

m=n,n' 



j^xx 



fc,m 



(18) 



The matrix on the RHS is the transpose of that in Eq. 
(11). Thus the two systems clearly have the same eigen- 



values, however, we prefer to use the w 



notation. 



together with its underlying creation and annihilation op- 
erators, because the overall wavefunction orthogonaliza- 
tions are determined by related commutators."'^^ A wave- 
function can be considered as a column vector of these 
coefficients: 



V'fe,2 



\ V'fe.N / 



(19) 



where some labeling from 1 to TV has been given to the 
N sites, each of which has a local two-component wave- 
function holding the S^ (w^ ) and S^ {w^} fluctuations: 



V'fe.n = 



W 



k.n 



(20) 



The dynamics of these individual wavefunctions can be 
written in terms of the matrices appearing in Eq. (18), 
which we denote here for convenience as A^nm: 



(21) 



m=n,n' 



Then we can see that there exists an operator A4 com- 
posed of the individual Ainm matrices, with real ele- 
ments, that evolves the entire wavefunction: 



itok'^k = M'^k- 



(22) 



This is the linear problem to be solved. 

We have previously solved this non-Hermitian prob- 
lem, Eq. (22), through numerical diagonalization, us- 
ing the EISPACK routine RG() for diagonalizing real- 
general matrices (no symmetry assumed). In this way 
the full spectrum could be found, however, the matrix 
M. is sparse but of size 2A'' x 2A'', where N is the number 
of lattice sites. A circular square lattice system with ra- 
dius R = 20a has N = 1264. This leads to a matrix with 
6390784 elements, which when stored in double precision 
(8 bytes per number), requires 49 MB RAM. In prac- 
tice, it is also necessary to store another array of equal 
size that contains the calculated eigenmodes, thus the re- 
quired storage for a calculation at one value of A is twice 
this number. This gets incremented by another 49 MB 
if one wants to compare eigenmodes at nearby values of 
A for tracking their changes with anisotropy. For these 
reasons, on a machine with 128 MB RAM, the upper 
limit for this calculation was system radius R =19-20, 
depending on the desired results. 

Next, we consider how the eigenmodes can be found 
through squaring the operator of the RHS of Eq. (22), 
combined with introducing a fictitious time evolution, re- 
sulting in considerable memory savings. 
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IV. RELAXATION METHOD 

A. Time Evolution 

A mode that solves Eq. (22) is a fundamental so- 
lution of 



dt 



(23) 



with time dependence e"^*^ , determined by the corre- 
sponding eigenvalue of Ai, iojk- Applying the M opera- 
tion again leads to an equivalent equation, 



(24) 



Both of these equations have fundamental oscillatory so- 
lutions varying as '^k{t) = 'l'fe(0)e*'^'=*, i.e., the usual 
spinwaves. The advantage is that operator has 
purely real, negative eigenvalues, in contrast to the purely 
imaginary cngcmvalucs (both positive and negative imag- 
inary parts) of operator Ai. 

Now consider that we invent a fictitious time evolution, 
by rather arbitrarily changing the LHS of Eq. (24) to a 
first time derivative: 



dt 



(25) 



The fundamental solutions here are determined by expo- 
nentials involving the eigenvalues of A4^, which are — w^. 
Thus, the general solution from an arbitrary initial state 
is 



*(t) = ^[c+*++c,-*,- 



-ojtt 



(26) 



where the coefficients Cf. arc determined by the initial 
state, and the are normalized eigenfunctions of M, 
with respective eigenvalues UtOk, corresponding to the 
creation and annihilation operators for mode k. (Of 
course, have the same eigenvalues, —lv^, with respect 
to the operator A4^.) All modes decay away in time, leav- 
ing the lowest frequency modes present in the initial state 
to dominate at large time. In practice, the time evolu- 
tion will be started from a randomly chosen initial state 
(i.e., randomly assigned values of the set of {w^,w^} 
coefficients), to ensure a nonzero overlap with the low- 
est frequency mode. The total wavcfunction ^'(t) will 
be re-normalized to unity periodically to avoid numeri- 
cal underflow, by enforcing the overlap appropriate for a 
non-Hermitian operator described below (Sec. IV B). 



(*|*) = 1. 



(27) 



We have used a second order Runge-Kutta time integra- 
tion scheme to solve Eq. (24) , and found that a time step 
of 0.02/ J S is adequate to provide stability and insure 
convergence to an eigenstate. 



At large time, the system will be relaxed to a linear 
combination of creation and annihilation wavefunctions: 



r.+ V]/ + 

-kn * fen 



(28) 



where ko labels the lowest mode that was present in the 
intial state. The frequency for this combination is evalu- 
ated in the usual way, 



(29) 



The creation and annihilation wavefunctions can be sep- 
arated by an application of operator M: 

M^ = iuJkoi4A-cg-^kJ, (30) 

Combining the two equations (28) and (30) leads to 



(31) 



This finally will be followed by a normalization appropri- 
ate to the non-Hermitian operator M, to be described 
below [Eq. (33)]. 



B. Orthogonality of Modes and Gram-Schmidt 
Process 

To obtain higher modes, it is necessary to evolve for- 
ward in time together with a Gram-Schmidt orthogonal- 
ization to any previously found lower frequency modes, 
denoted here by index i. We have found that the neces- 
sary way to do this is by making the current wavcfunc- 
tion ^'(t) orthogonal to both and rather than to 
a somewhat arbitrary combination such as in Eq. (28). 
Equivalently, both constants and cj must be enforced 
to be zero in the current wavefunction, for all modes £ 
already determined. To obtain these two conditions, it 
is necessary to use the overlap appropriate to the non- 
Hermitian operator A4, since Vl*^ and 'ij are distinct 
eigenfunctions of Ai corresponding to its distinct eigen- 
values ±iuJi. 

The overlap of two eigenmodes of M can be defined 
in one of two equivalent ways. The first is to realize 
that an annihilation {Bj) and creation {bI) operator, 
with associated wavefunctions '^J and ^'l, must have a 
commutator, 



(32) 



Prom the definition, Eq. (15), with the help of the basic 
spin commutators [5'n,5'^/] = iSn^n'S^ « iSn^n'S, this is 
equivalent to: 



[Bj,Bl] 



2 * 



(33) 



5 



In general, this is the definition of the overlap used be- 
tween two arbitrary states j and k. 

The second way to define the overlap, equivalent to 
this, is to make use of the left and right eigenvectors 
of M. The solutions of Eq. (22), where M acts to 
the right, are the right eigenvectors. For each eigenvalue 
iojk, there also exists a corresponding left eigenvector, fc\E' 
satisfying 



(34) 



Each creation operator as well as each annihilation op- 
erator has a distinct left eigenvector. For a Hermitian 
problem, the left eigenvectors are obtained by the famil- 
iar relationship, f.'i/ = 'i/j., and then overlaps are in the 

usual form {^j\^k) = j^'^fc = 

For this problem, although Ai is not Hermitian, there 
exists a real matrix A= —A~^, that transforms A4 to 



where A is diagonal in 2 x 2 submatrices:^^ 



-1 

1 



(35) 





/ a 





. 











a 


. 






A = 








a . 




a 




I '. 






'.'.) 





(36) 



Then using Eq. (35) in Eq. (22), rearranging and com- 
paring with Eq. (34) determines a relationship between 
left and right eigenvectors, 



fc* = i (^*fc)'^ 



(37) 



where the f operation is the usual Hermitian conjugate 
(transpose of complex conjugate), and the factor of i has 

been added for consistency with Eq. (33) . Then the over- 
lap between two states and "ifk is defined by the scalar 
product of a left with a right eigenvector: 

(*,|*fe) = = I (A^,)^ *fe = z^-j^^^fc. (38) 

which is seen to give exactly the expression in (33). 
The application of Gram-Schmidt orthogonalization to 

remove the previously found modes (labeled by £) pro- 
ceeds according to replacing the current "^{t) by 

*(t) - *(t) - ^ [(*+|*(t)) *+ + <vi/7|*(t)> *7] . 

e 

(39) 

Combined with evolution forward in time and periodic 

renormalization according to Eq. (27), the wavefunction 
will evolve to the next low-frequency mode, of the form 
Eq. (28). Continuing to iterate the above procedures, a 
small set of the lowest frequency modes (about 5-10) can 
be determined with double precision accuracy in reason- 
able time for systems of size 100 x 100. 



C. Self-Consistent Gauss-Seidel Scheme 

The eigenmodes of Eq. (22) alternatively can be found 

by a Gauss-Seidel type of iteration, not applied to that 
equation, but rather, to the equation with the squared 
operator, 



n = M^ 



(40a) 



(40b) 



This Gauss-Seidel scheme to be described is about twice 

as fast (in total cpu time) as the time evolution using 
second order Runge Kutta described above. It is based 
on using the matrix elements of H, and in particular, 
one needs to single out the matrix elements that couple 
a site to itself (Wn.n) from the intersite couplings (Wn). 
We write this equation as a set of equations with the 
form, 



22 



W. 



n" /3=1,2 



/3 



(41a) 
(41b) 

(41c) 



where the sites n" include all first and second nearest 
neighbors of site n, but not the site n. The second near- 
est neighbor interations result directly from squaring the 
nearest neighbor operator, M.. These matrix elements of 
H are sums of quadratic products, 



E 



(42) 



m=n,n' ,n" 



where the summation is over any sites m that are either 
first neighbors of sites n and n", or, these sites them- 
selves. Due to the structure of the matrix elements of 
A^, the on-site, ofF-diagonal elements 7in^„, T^n^-n, are 
both zero, while the corresponding diagonal elements are 
equal: = 

Now in the spirit of Gauss-Seidel iteration, we can try 
to "solve" for the new values, i/^n', with some certain 
freedom of choice. That is, the terms involving can 
be taken as explicit or implicit, depending on whether we 
consider them to involve the wavefunction at the previous 
step or the wavefunction at the current step being solved. 
We have found for this problem, that the iteration is more 
stable if the terms are explicit, and the equation is re- 
arranged for each component a = 1,2 as 



-1 



(43) 



For each step of the iteration, the frequency lo vaxist be 
evaluated according to an expectation value of the form 
of Eq. (29). Also, in the interest of numerical efficiency, 
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wc actually store the matrix elements of A4 and only the 
diagonal elements of 7Y = A4^ . Then it is useful to realize 
that an equivalent expression for Wn is in terms of the 
entire H operation with the diagonal part subtracted out: 

ws = [mt-n^>^] (44) 

This means Eq. (43) can be written as 

<' = <-^K< + (W*)a- (45) 

'*'n,n 

The stability of this process is further enhanced if we mix 
a fraction / < 1 of the updated configuration [Eq. 
(45)] with a fraction (1 — /) < 1 of the old configuration 
^. This then finally defines the iteration procedure: 

wZ' = wZ-:^[u;'wZ + in^t]. (46) 

The equation is to be iterated, together with Eq. (29) for 
the frequency and the usual normalization, Eq. (27), until 
a self-consistent convergence is achieved. We should note 
also that the application of Eq. (46) can be performed ei- 
ther synchronously or asynchronously. In a synchronous 
step, the old values «;„ arc updated into new ones w^' 
simultaneously, or in parallel. In an asynchronous step, 
which is the usual Gauss-Seidel scheme, once a site is up- 
dated into the new value w^'i that new value is recycled 
immediately into the right hand side of Eq. (29), and 
this new value then modifies the results for its neighbor- 
ing sites. Either way, the value of lo'^ is corrected only 
after the entire lattice has been updated. Both ways work 
well, although the asynchronous step generally produces 
faster convergence. Just as in the time evolution scheme, 
a random initial configuration of ^' is used for these iter- 
ations. In our actual calculations, we have used fractions 
/ from 0.6 to as high as 1.9 (an over-relaxation). Values 
of / below 1.0 provide reasonably fast convergence with 
reliable stability. Fractions / that arc closer to zero result 
in slower convergence, while fractions closer to one give 
faster convergence at the expense of greater instability, 
where the instability causes convergence to the highest 
frequency eigenmodc rather than the desired lowest fre- 
quency mode. 

The above procedure will tend to the lowest frequency 
mode. As for the time evolution scheme, to get higher 
modes it is only necessary to combine the Gauss-Seidel 
iteration with the Gram-Schmidt orthogonalization pro- 
cedure described in Sec. IV B. 



D. BoundEiry Conditions 

To completely define the eigenvalue problem, bound- 
ary conditions must be specified. The simplest choices 
arc Dirichlet and free boimdary conditions. For Dirichlet 
boundary conditions, we imagine placing fictitious spin 
vectors Sn,,, fixed to the static vortex directions, at sites 



nb just outside the system. The sites just inside the sys- 
tem interact with these, leading to on-site terms in the 
matrix M due to the boundary, see Eq. (13b), 

m2 = -Mil 

= JSY, [P>°n, cos(0O - 00 J + AmOmO J . (47) 
"b 

These arc the only effects of the boundary. 

For free boundary conditions, the system is open at 
the boundary, and no fictitious spins are needed. There- 
fore, the above terms of Eq. (47) arc absent. This is the 
only difference in the matrix between Dirichlet and free 
boundary conditions. 

However, in the absence of an in-plane magnetic field, 
free boundary conditions lead to a zero frequency mode, 
due to the in-plane rotational invariance of the model. 
This mode is unusual in that it cannot be made orthogo- 
nal to the other modes by the overlap integrals discussed 
above. Therefore, it needs to be removed from the evolv- 
ing mode solution by a different tactic. 

Consider its structure for A < Ac, i.e., in-plane vortices, 
where = 0, p° = 1 on all sites. Eqns. (18) become 
fairly simple, 

JS cos{<l,l-<Pl){wi-w'J=iwkwi,, (48a) 

m— n,n' 

JS Yl {cosicPl - cl>M - Xw^} =iu;kwl (48b) 

m— n.n' 

For the zero frequency mode, there is the solution, 

where <fo is an arbitrary constant. Obviously this mode 

is not normalizable by (33), and represents the fact that 
a uniform in-plane rotation costs no energy. 

For A > Ac, i.e., out-of-plane vortices, there is no great 
simplification of (18), however, a short calculation shows 
that the zero frequency mode is 

*'-=(<:)^(...™<). M 

where ipo again is an arbitrary constant, and 0° is the 
static out-of-plane structure of the vortex. This wave- 
function is slightly diminished in amplitude near the vor- 
tex core, where the spins point out of the easy plane. 
One can think of the component as a "coordinate" 
{w^ = <PnCOS^° w in-plane angle), which has the spec- 
ified value, and W'' as the corresponding conjugate mo- 
mentum (part of due to the spinwave), which has been 
set to zero. 

In the evolution to find other modes, any component of 
the current solution proportional to V'o must be removed. 
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For A < Ac, it is fairly obvious that we only need to en- 
force the constraints. = 0, (w^) = 0, where ( ) 
indicates the spatial average. These are essentially the 
conditions that the total conjugate momentum is zero, 
and that modes have no uniform in-plane drift compo- 
nent. 

For A > Ac, physically, the higher modes must also 
have a total conjugate momentum equal to zero, and, 
clearly, must satisfy {(f)n) = 0. This leads to the general 
constraints we applied for free boundary conditions, 

(55^) = K cos 0=0, (51a) 

{fn) = {w';,/cosel) = 0, (51b) 

where we used 3^ = 8^" + SS^ = sin 6*° -I- cos 0^^ to 
define the dynamical conjugate momentum, 6S^. 

E. Memory and CPU Considerations 

These relaxation methods use considerably less mem- 
ory than the diagonalization scheme, provided that we 
arc satisfied to obtain only a few of the lowest frequency 
modes, which is the case here. In order to obtain good 
stability for systems up to radius R = 100a, it is neces- 
sary to use double precision. 

If we want to calculate the lowest K modes of a sys- 
tem N sites, then we need to save 4N complex numbers 
for each state to be calculated (creation and annihilation 
operators), or 8NK doubles for all states. The matrix el- 
ements of M are real and occupy (4 x 4 + 2) x TV = 18A^ 
memory locations, where the factor accounts for not stor- 
ing the elements that are zero. For either type of re- 
laxation scheme, it is convenient to have two complex 
arrays each with 2N elements to hold M'^ and M'^'if 
{8N doubles). For the Runge-Kutta time evolution, an- 
other two complex arrays of size 27V are needed. For the 
Gauss-Seidel iteration, on the other hand, the real terms 
^733- are saved, requiring N memory locations. For the 
asynchronous Gauss-Scidel, cpu time is reduced by sav- 
ing Hn,n", which requires 4 x 12 x iV doubles (there is 
a total of 12 nearest and next-nearest neighbors on the 
square lattice). Thus the total program size required for 
a calculation at a single value of anisotropy A scales as 
(34 + 8K)N for Runge-Kutta and as (27 + 8K)N (syn- 
chronous) or as (75 -|- 8K)N (asynchronous) for Gauss- 
Seidel. 

For a circular system of radius R = 20a, with N = 
1264, and choosing K = 15, these come to 1.52 MB 
for Runge-Kutta and 1.45 MB or 1.93 MB for Gauss- 
Seidel, considerably smaller than the approximately 100 
MB necessary by diagonalization. The required memory 
increases if we want to track the changes in the modes 
with changing anisotropy, which requires storage of the 
set of modes at two different values of A. The primary 
drawback of these methods is that the cpu time may be 



quite long, because the solution method is essentially dif- 
fusive. We find that the cpu time increases as R*, where 
one factor of R^ is the increase in system area, and the 
other factor of R^ is due to the diffusion time increasing, 
tdiffuse oc R^ . Often we are interested to calculate for a 
set of 20 to 40 values of A between zero and one; the large 
cpu time has prevented us from going beyond system ra- 
dius R = 100a, where each mode may take on the order 
of a day on typical 100 MHz processors. 

It is also clear that the rate at which each mode con- 
verges diminishes with increasing system size, because 
the frequency differences between modes decreases as 
1/R. Furthermore, it is especially difficult to converge 
to any mode that is barely split from the next higher 
mode. However, the doubly degenerate modes present in 
this magnetic spinwave problem do not cause difficulties. 
Any arbitrary linear combination of such a pair will solve 
Eq. (22), and then the Gram-Schmidt process during the 
relaxation to the next higher mode will always produce 
another equal frequency mode, orthogonal to the first, 
spanning the degenerate subspace. 

V. RELATION TO CONTINUUM THEORIES 

For a full theoretical description of vortex-spinwave in- 
teractions we need to mention how to compare our results 
for the wavefunctions with the spin fields determined in 
continuum limit theories. This is straightforward for the 
FM model, for which we will only make a few comments. 
For the AFM model, however, we are usually analyzing 
the staggered magnetization i. Below we describe how i 
is related to the fields. 

A. FM Continuum Coordinates 

For the calculations on a lattice described so far, it 
was most convenient to use the Cartesian spin compo- 
nents. To make some comparison to theory wc used pla- 
nar spherical coordinates [Eq. (6)] , with ^pianar measured 
from the ajy-plane, such that ^pianar is conveniently zero 
in the ground state. Of course, in much of the literature it 
is common to use polar spherical coordinates, with f?poiar 
measured from the positive 2:-axis. Clearly, these are re- 
lated by 0poiar = 7r/2 — ^planar- Then the relation that 
replaces (10), between our local Cartesian components 
and the angular fluctuations, using polar spherical coor- 
dinates, is 

= 5(p„ sin ^° , Sl = -5i?„. (52) 

The reason wc mention this is that the quantity /i — 
(p{r) sin0°(r) appears naturally in the continuum theory 
for FM vortex-spinwave scattering,'' and it is equivalent 
as well, via Eq. (17), to our field. Similarly, our 
field represents the i?-fluctuations discussed in continuum 
theory. 
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B. AFM Continuum Coordinates 

For the AFM scattering problem, there have been two 
primary ways to approach the continuum description, 
where we consider only lattices composed of two sub- 
lattices. 

For example, in Rcfs. 10, 12, the theory is developed 
directly in terms of £, defined by the difference of spins 
on the two sublattices, 

i=^^LZl^=S^ + ^^±^i^ + ^^zJky^, (53) 



25 



25 



25 



where the local axes are related by Zn' = —Zn, in' = 
—Xn, yn' = j/n, and n, n' are neighboring sites, on the 
two different sublattices. Assuming £ has perturbations 
and is described by planar spherical angles (pi = (p'j + ipg, 
0e = 9^ + 'di, and using Eq. (17), we then obtain the 
equivalences. 



25 



sj-sj, 

25 



25 



-wt, 



ifccosOl, (54a) 



25 



(54b) 



Here the first cqiiation gives the in-plane fiuctuations 
while the second gives the out-of-plane tilting fiuctua- 
tions. The mapping to the fields is correct up to an 
arbitrary phase. If polar spherical coordinates are used, 
then change cos^^° smO^^ and -di Once again, 

it is found that is equivalent to the variable /i£ intro- 
duced by Ivanov et al,^^ that appears together with 
in symmetrical Schrodinger-like equations. 

Alternatively, in other vortex-magnon scattering 
articles, the spins on the two sublattices have been 
represented in polar spherical coordinates introduced by 
Mikeska,^^ using four angular fields, as 



sin(e ± 9) cos($ ± (/)) 
±5 I sin(e ± 9) sin($ ± ^) 
cos(e ± 9) 



(55) 



where for simplicity we suppressed the site indices, and 
upper/lower signs refer to the two sublattices. The 
"small" angles 9 and <p arc slave to the "large" angles Q 
and via 9 = sin 6/8 J5, and = e/(8J5sine), 
showing that both are truly small for low-frequency 
modes. Then assuming that the large angles have small 
perturbations away from a vortex configuration 8°,$°, 
O = + 1?, <i> = + a short calculation show 
that the staggered magnetization will have components, 
= (fisinQ^, = We see that these are equivalent 
to the perturbation variables in Eq. (54) , once the change 
from polar to planar coordinates is taken into account. 



VI. APPLICATION: VORTEX-SPINWAVE 
SPECTRA AND S-MATRIX 

The methods described above were tested by compar- 
ing their results with numerical diagonalization results^ 
for small systems (up to radius R w 20a). At these small 
system sizes, it was even possible to get very precise re- 
sults (more than 6 digits for eigenvalues) for as many as 
the lowest 50 modes. Typically, for A / 0, the mode fre- 
quencies for the AFM model are higher than those for 
the FM model, which means that the calculations tend 
to converge faster for the AFM model. 

As an example of the utility of these methods, we con- 
sider the spinwave spectrum (lowest 30 modes) for a sin- 
gle vortex in a circular system, with Dirichlet boundary 
conditions, and use the results to calculate the scattering 
S-matrix. We used a square lattice of spins. For most 
of the data presented, the synchronous self-consistent 
Gauss-Seidel scheme with fraction / = 0.7 was used. 
The FM model for A very close to the critical value is an 
exception; there we used the asynchronous sheme, with 
/ = 0.9, which has greater stability at the expense of 
greater CPU time. 

Typical mode spectra for some small systems were 
given in Ref. 2, as functions of the anisotropy param- 
eter A. In the present work we can consider system radii 
as large as i? = 100a, however, the CPU time to get the 
complete dependence versus A at such a large radius is 
prohibitive, and for most calculations we used R < 50a. 
We concentrate on the spectra at a few isolated values 
of A < Ac, corresponding to in-plane vortices, and focus 
more on the variations with system radius R. 



A. FM/AFM In-plane vortex S-Matrix 

Due to azimuthal symmetry, the modes generally are 
classified according to azimuthal quantum number m, 
corresponding to an expected e*"*^ dependence, where 
a point in the system has polar coordinates (r, x) , mea- 
sured from the center. Similarly, a principal quantum 
number n can be defined, that is the number of nodes 
in a radial direction in the wavefunction. There are only 
minor differences in calculating a scattering amplitude 
Pm{k) [defined below, see Eq. (56)] for the FM model 
compared to the AFM model. For FM out-of-plane vor- 
tices, the rotational symmetry around the 5^ axis is bro- 
ken by the nonzero topological charge and nonzero 
spin components, and generally, m and —m modes are 
not degenerate. For the in-plane case we consider here, 
however, the topological charge is zero, rotational sym- 
metry holds, and ±m modes are degenerate. In the con- 
tinuum limit, the zero-field AFM model has complete 
rotational and inversion symmetry around the 5^ axis, 
with no topological charge, so that the modes ±m are 
degenerate at any value of A. For both FM and AFM 
in-plane vortices, the lack of 5^ spin components insures 
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degeneracy of ±m modes. The degeneracy is partially 
broken on a lattice;^ for even m, the eigenfunctions gen- 
erally are mixtures of +m and —m components. 

We can analyze the scattering S-matrix both for 
FM and AFM models as follows. In the continuum 
description/^ the lowest modes, which belong to the 
acoustic spinwave branch, have the asymptotic form, 



w 



'(r-,x)«0, 



(56a) 



u.2(r,x) - [J™(fcr) +p„(fc)y™(fcr)]e^"^ (56b) 

where Jm and Ym are Bessel and Neumann functions, and 
k is the wavevector as determined from the frequency 
of the mode. For the AFM model, the free spinwave 
dispersion relation is 



Wk = 4J5v/(lT7k)(l±A7k), 



(57) 



where upper/lower signs are for the acoustic/optical 
branch, and 



7k = -^{cosk^c + cosky). 



(58) 



For the FM model, we have only an acoustic branch with 
dispersion relation 



= 4J^v/(l-7k)(l-A7k). 



(59) 



Note that the AFM acoustic branch dispersion can be 
obtained from the FM dispersion simply by the change 

A —A. After the relaxation scheme is used to obtain a 
mode frequency lu, then we invert the dispersion relation 
to get the corresponding wavevector needed in the fitting 
expression (56b). For the FM dispersion, assuming k = 
{k,0), we have k = cos~^(27k — 1), where 

7k = (2A)-i {(1 + A) - v'(l + A)2-4A[l-(a;/4J5)2]} . 

(60) 

The same formula applies to the AFM acoustic branch 
when the change A — » —A is made. 

In Eq. (56), the coefficient Pm{k) is a measure of the 
scattering; in the absence of scattering it vanishes, and 
the free spinwave is described by the Jm{kr) function 
alone. Alternatively, the wavefunction can be ex- 
pressed using the S-matrix scattering function, Sm{k) = 
g2iA„(fe)^ where Am{k) is the phase shift, in terms of 
incoming and outgoing waves, as 



w 



2'Kkr I 



-i{kr+4'm) 



+ Sm{k)e 



i{kr-\-<prn) 



}, (61) 



where the phase angle is (f)m = f (2m + 1). In the ab- 
sence of scattering, Sm{k) ^ 1, or Am(fc) 0, and the 
asymptotic form of Jm{kr) results. Expressions (56) and 
(61) are equivalent up to a constant, provided 



Sm{k) — 



ipm{k) 



1 + ipm{k) ' 
Then the phase shift is related to Pm{k) by 
tanA„(A:) = -pmik). 



(62) 



(63) 



The scattering amplitude pm{k) can be found by mak- 
ing a least squares fitting of a wavefunction to the above 
asymptotic form (56b). However, in the case to > 0, 
when TO and —to are degenerate, we actually fit to the 
expression, 



w\r,x) = [aMkr) + WY^{kr)] e^™^ 
+ [a2Jm{kr) + h2Ym{kr)] e-'^^, 



(64) 



where ai, a2, hi, and 62 are complex fitting constants. 
Then, Pm{k) = bi/ai and p^m{k) = 62/^2 arc extracted, 
and should be expected to be comparable. For in-plane 
vortices, whose static structure has no A-dependence, we 
considered the asymptotic region to be r > 8a. 



B. The Scattering Results: In-Plane Vortices 

We considered in-plane FM and AFM vortices for 
A = 0.0, 0.5, 0.7. The static in-plane vortex spin structure 
is confined to the xy plane, resulting in a singular point 
at the vortex core, and making continuum theory diffi- 
cult. For any A < Ac, the static structure of in-plane FM 
vortices is the same as the in-plane AFM vortex structure 
on one sublattice. However, the dynamics are different, 
in particular, the FM mode frequencies tend to be much 
lower than those for the AFM, although the appearances 
of their wavefunctions are very similar. The only excep- 
tion to this is at A = 0, where the FM and AFM models 
have the same mode frequencies. Thus it is interesting 
to compare the FM/AFM scattering properties and also 
investigate the anisotropy dependence. 

Wavefunctions of some of the lowest modes were pre- 
sented in Refs. 12, 2, and more complete wavefunc- 
tion and other data, from the new methods described 
here, can be found at http://www.phys.ksu.edu/"' 
wysin/vortexmodes/ . From those diagrams it is easy to 
identify the principal and azimuthal quantum numbers n 
and m for each mode. In Fig. 1 we show the dependence 
of some of the lowest eigenfrequencies on system radius 
R, for A = (identical results for FM and AFM here). 
Differences appear in these spectra as A approaches the 
critical value, as seen in Fig. 2, where the spectra are 
shown for A = 0.70, which is just below Ac ~ 0.7034 . 
Additional data for A = 0.5 are presented at the www 
page listed above. With the exception of the local mode 
in the AFM, most modes have frequencies diminishing 
asymptotically as 1/R, as expected in a continuum the- 
ory. The local mode in the AFM model has a frequency 
independent of the system size, but dependent on A, and 
can be thought of as a state taken out of the optical spin- 
wave branch. The lowest mode in the FM model is only 
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quasi-local: its wavcfunction is local in the out-of-plane 
fluctuations but extended in the in-plane fluctuations, 
and its frequency decreases slowly with increasing sys- 
tem size. 

We calculated the spectrum of the lowest 30 modes 
for different system sizes R, ranging from R = 13a to 
R = 100a. Inverting the ujk relationship (59) as described 
above [via Eq. (60)], and by using these different sys- 
tems sizes, we could observe a particular mode labeled 
by (n, to) at different k. In this way we obtained the 
fc-dependence of the scattering. 

1. Scattering by FM Vortices 

Some results obtained for the scattering of spinwaves 
by FM vortices are shown in Fig. 3 and Fig. 4, where the 
scattering amplitude Pm{k) is displayed for A = 0.0,0.7. 
The scattering tends to be much weaker as m increases 
above 2. As fc — > 0, pm{k) tends to zero, a physically rea- 
sonable limit. An interesting result, especially as A ^ Ac, 
is that Pm{k) has points where it is singular and changes 
sign. Indeed, the to = scattering curve at A = 0.7 ex- 
hibits two singularities, one near ka = 0.035 and another 
near fca = 0.65. Note that the scattering results obtained 
for A = 0.0,0.5, probe down to ka = 0.024; down to this 
point they do not display the lower singularity. 

Using Eq. (63), we see that these singular points 
in Pm{k) correspond to the phascshift A„i(fc) passing 
through the value 7r/2 near ka = 0.035 and near ka — 
0.65. In Figs. 5 and 6 we show the resulting phaseshifts 
^mik) for A = 0.0,0.7; the phaseshifts are smooth func- 
tions oik. It is impressive that the to = phase shift for 
A ~ 0.7 rapidly surpasses 7r/2, reaches a maximum and 
then falls back through it/ 2. Additional data for other 
parameters is found at the www page. Note that as long 
as Pm{k) <C 1 there is little difference between it and the 
phase shift, except for the reverse sign. 

Costa et al.^ calculated the to. = and to = 1 phase 
shifts for the FA4 model, only for A = 0, in the contin- 
uum limit via a Born approximation. Our results here are 
somewhat in contrast to theirs; from Fig. 3 and Eq. (63) 
one sees that we obtain the opposite sign, (k) > for 
to = 0, 1, 2, and over the range, < fca < 0.5, there are 
no singularities in Pm{k). Thus these phase shifts do not 
reach tt/2 in this interval. For to = 0, the Costa et al. 
result is similar to ours (except for the opposite sign), 
whereas for to, = 1 their result for the phase shift Ai(fc) 
has reached — tt at ka « 0.3. Note, however, that the ana- 
lytic calculations require the (somewhat arbitrary) choice 
of a short distance CTitofF on the integrals appearing in 
the Born approximation, at a distance on the order of 
a lattice constant (a/2 was used in Ref. 8). It is clear 
that a different cutoff could lead to drastically modified 
results, as was shown, for example, in a calculation^^ of 
the frequencies of the local mode of AFM in-plane vor- 
tices. The choice of how to smoothly convert the singular 



field in the core of the vortex into a smooth continuum 
field is not clear, and always makes it difficult to write 
out an unambiguous continuum theory for the in-plane 
vortices. Therefore it may not be so fruitful to make any 
detailed comparison of the discrete lattice results here to 
continuum results for in-plane vortices, unless no cutoff is 
used. (This is not a problem for the out-of-plane vortices, 
which are not singular at their core.) 

2. Scattering by AFM Vortices 

At A = 0, there is no difference in the scattering for 
AFM vortices compared to FM vortices; Figs. 3 and 5 ap- 
ply. (Only the acoustic AFM branch is being considered; 
the optical branch is too high in the spectrum for these 
relaxation techniques to be useful.) Similarly, results for 
the scattering by AFM vortices at A = 0.7 are shown in 
Fig. 7. Results obtained at A = 0.5 are very similar to 
these, and shown on the web page indicated above. For 
the to = 1 scattering, there is a singular point somewhere 
around ka w 0.55; thus Ai(fc) passes through 7r/2 at this 
point, as seen in Fig. 8. Although the AFM frequencies 
are higher than the corresponding frequencies for the FM 
model, the data obtained probe to similar low values of 
ka « 0.024, which is essentially determined by the largest 
system size used, R = 100a; no singularity for to = was 
found down to this wavevector limit. Therefore we can 
say that the m = scattering by AFM in-plane vortices, 
just below Ac, is considerably weaker than that for FM 
in-plane vortices. On the other hand, the m — 1 scat- 
tering in this case is much stronger by the AFM in-plane 
vortices than by the FM in-plane vortices. 

Continuum results for the AFM acoustic branch spin- 
wave scattering by in-plane vortices are intriguing. 
Pereira et al}^ have calculated the scattering phase shifts 
for the optical branch, which excite o\it-of-plane spin fluc- 
tuations (our field), even without the need for a Born 
approximation. However, in the leading order continuum 
theory, the potential due to the vortex appears only in 
the perturbation equation for the out-of-plane fluctua- 
tions (angle dg); it is absent from the equation for the 
fluctuations of the in-plane components (angle (pi or our 
field). Indeed, the dynamical equation for ipg is sim- 
ply a wave equation for free acoustic spinwaves. It means 
that the usual continuum limit predicts no scattering of 
acoustic spinwaves (which excite in-plane fluctuations) 
by in-plane AFM vortices. We might comment, however, 
that for out-of-plane vortices, the potential due to the 
vortex appears in both the equations, one for the 
other for the quantity pe = (pf cos 6'^ [Eq. (54)]. So the 
scattering from in-plane vortices is a special case. 

Of course, with our numerical calculations on the lat- 
tice for AFM in-plane vortices, we did obtain some rela- 
tively weak scattering for to = 0, 2, and rather stronger 
scattering for m = 1. It is possible that this unexpected 
scattering is caused primarily by the vortex core, where 
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the static spin directions rotate by 90° for nearest neigh- 
bor sites. This gives large gradients in the core region, 
which cannot be described in the usual lowest order con- 
tinuum limits. Some theoretical analysis of this kind of 
discrete effect, which always occurs around the vortex 
core, is needed. 

VII. VORTEX MASS 

As another application, information about certain nor- 
mal modes can be used to estimate a vortex mass. The 
results here apply only to in-plane vortices. The modes 
with (n = 0,m = 1), which are doubly degenerate for 
in-plane vortices, are most strongly associated with the 
translational motion of the vortex center. If only these 
modes are added to the original static vortex structure, 
the resulting configuration produces either a circular or 
linear simple harmonic oscillatory motion of the vortex 
center. -"^'-"^^ These motions are the same as those predicted 
from a simple Newtonian dynamical equation of motion 
for the vortex center, involving a force F and effective 
mass M, i.e., 

F = MV, (65) 

where V is the velocity of the vortex center X{t). The 
translational modes can be considered to be driven by 
the force caused by the system boundary and the lat- 
tice itself, which can be evaluated independently of the 
spinwave spectrum. In the simplest approximation, we 
assume a linear restoring force, F = —KX, where the 
spring constant K has been found by a relaxation scheme 
applied to a static vortex in Ref. 14. For in-plane vor- 
tices, it was found that K w AJS'^/a'^, independent of 
the value of A. Then, for this simple harmonic motion, 
the mass is found as 

M = K/ujI^, (66) 

where Wo,! is the translation mode frequency. The re- 
sults of this calculation arc shown in Fig. 9, for both 
the FM and AFM models at various anisotropics. The 
fact that the translational mode frequencies diminish as 
the reciprocal system size is reflected there; the effec- 
tive mass increases as the square of the system radius R. 
The masses are identical for FM and AFM vortices at 
A = 0, where their spinwave spectra are identical. The 
FM vortex mass increases with increasing A, while the 
AFM vortex mass decreases. Thus we have the result 
that in some real dynamical sense AFM vortices will be 
easier to move than FM vortices. The fact that the mass 
depends on the system size should not be unexpected, 
because a vortex is not a localized object. Its spin con- 
figuration extends to the limits of the system, causing 
its inertia against motion inside the system to depend on 
the system size. The mass defined here really is more a 
property of the entire system. 



VIII. CONCLUSIONS 

We have discussed some relaxation methods for deter- 
mining the lower part of the eigenspectrum for magnetic 
models in the presence of an individual magnetic nonlin- 
ear excitation, such as a vortex. The methods are based 
on the fact that the effective time evolution of the squared 
Hamiltonian is diffusive, and therefore evolution from an 
arbitrary intial state will lead to the lowest mode present. 
By combining with Gram-Schmidt orthogonalization, a 
reasonably large set of the lowest modes, including their 
eigenfunctions and eigenfrcquencics, can be determined. 
Special considerations also were described for applying 
free boundary conditions, in contrast to Dirichlet bound- 
ary conditions; for the former the zero- frequency mode 
must be eliminated from the evolving solution. 

As an application of the schemes, the vortex-spinwave 
scattering and S-matrix were determined numerically for 
both FM and AFM in-plane vortices on square lattice cir- 
cular systems. The results obtained exhibit some inter- 
esting singular points in the scattering amplitude pm{k), 
especially as A approaches the critical value from below. 
A singularity in the rn = 1 scattering amplitude by FM 
or AFM vortices (A = 0.0) occurs near ka « 0.65. A sin- 
gularity in the m = scattering amplitude for the FM 
at A = 0.7 occurs at very low wavevector (ka « 0.035) 
or over long distances, and it may be interesting to ask 
whether it could be described in a continuum theory. 
These singular points occur where the phase shift passes 
smoothly through values like ±7r/2. 

At A = 0, the scattering results for FM and AFM vor- 
tices are identical, as a result of identical mode frequen- 
cies. For the AFM model, although theoretically there 
should be no scattering of acoustic spinwaves by a vortex, 
we did find nonzero scattering amplitude, especially for 
m = 1, and even a singular point in this amplitude for 
A = 0.7 . This imcxpected strong scattering is probably 
due to the large gradients in the spin field around the 
vortex core, which cannot accurately be described in the 
usual continuum limit theory. 

The lowest vortex translational modes (lowest mode 
most closely coupled to motion of the vortex center) were 
used to estimate a vortex mass, for in-plane vortices, in a 
very simple phenomcnolgical theory. The theory is based 
on the idea that the vortex in a small system vibrates in 
its translational mode in response to the force on it due 
to the lattice and system boundary, with a frequency 
determined by the effective force constant and the vortex 
mass. Using the observed translational mode frequencies, 
this simple theory leads to a vortex mass increasing as 
the square of the system size, for both FM and AFM 
vortices. AFM vortices have lower effective mass and 
thus are expected to be easier to move. For the in-plane 
vortices, FM vortices at A ~ Ac have the highest mass, 
while AFM vortices at A ~ Ac have the lowest mass. 

With appropriate modifications to the Hamiltonian, 
etc., in general, the relaxation schemes discussed here 
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would apply to other kinds of stable nonuniformly mag- 
netized systems, such as, for example, small magnetic 
particles with a frozcn-in magnetic configuration. Then 
the methods might be very useful for determining any 
interesting response properties of such materials. 
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FIG. 1. The frequencies of some of the lowest modes for a 
single FM or AFM in-plane vortex, at A = 0.0, versus sys- 
tem radius R. The integers (n, m) label the principal and 
azimuthal quantum numbers. 
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FIG. 3. The lowest scattering amplitudes Pm(fc) for a FM 
or AFM in-plane vortex, at A = 0.0, for a) m = 0,1,2 and 
b) expanded view of the m = 1 amplitude, exhibiting the 
singularity. 
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FIG. 2. The frequencies of some of the lowest modes for a 
single a) FM or b) AFM in-plane vortex, at A = 0.7 < Ac, 
versus system radius R. 
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FIG. 4. The lowest scattering amplitudes pmik) for a FM 
vortex, at A = 0.7 < Ac, for a) m = 0, 1 and b) m = 2. 
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FIG. 5. The lowest scattering phase shifts Am{k) for a FM 
or AFM in-plane vortex, at A = 0.0, for m = 0, 1, 2. 
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FIG. 6. The lowest scattering phase shifts Am{k) for a FM 
in-plane vortex, at A = 0.7, for m = 0, 1, 2. 



FIG. 7. The lowest scattering amplitudes pmik) for an 
AFM vortex, at A = 0.7, for a) rn = 0,2 and b) m = 1, 
showing the singularity between ka = 0.5 and ka = 0.6 . 
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FIG. 8. The lowest scattering phase shifts Ajnik) for an 
AFM in-plane vortex, at A = 0.7, for m = 0, 1, 2. 
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